NPS55-90-20 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


LATENT  FACTOR  MODELS  AND  ANALYSES 
FOR  OPERATOR  RESPONSE  TIMES 


Donald  P.  Caver 
I.  G.  O'Muircneartaigh 


September  1990 


FedDocs 

D   208. 14/2 

NPS-55-90-20 


Approved  for  public  release;  distribution  is  unlimited. 

Prepared  for: 

Naval  Postgraduate  School 
ey,CA  93955 


NAVAL    POSTGRADUATE    SCHOOL 
MONTEREY,    CALIFORNIA 


Rear  Admiral  R.  W.  West,  Jr.  Harrison  Shull 

Superintendent  Provost 


This  report  was  prepared  in  conjunction  with  research  funded  under  the 
Naval  Postgraduate  School  Research  Council  Research  Program. 

This  report  was  prepared  by: 


iqiJj  M9Y 


SCHOOL 


■sification  of  this  page 





REPORT  DOCUMENTATION  PAGE 


security  Classification  UNCLASSIFIED 


1  b    Restrictive  Markings 


Classification  Authority 


fication/Downgrading  Schedule 


3       Distribution  Availability  of  Report 

Approved  for  public  release;  distribution  is  unlimited 


g  Organization  Report  Number(s)    NPS55-90-20 


5     Monitoring  Organization  Report  Number(s) 


Performing  Organization 

tgraduate  School 


6b  Office  Symbol 
(If  Applicable)    OR 


7a  Name  of  Monitoring  Organization 


(city,  state,  and  ZIP  code) 

CA  93943-5000 


7  b  Address  (city,  state,  and  ZIP  code) 


Funding/Sponsoring  Organization 

Postgraduate  School 


8b  Office  Symbol 
(7/  Applicable) 


9     Procurement  Instrument  Identification  Number 

O&MN,  Direct  Funding 


(city,  state,  and  zip  code)  Monterey,  California 
;y,  CA    93943 


1  0    Source  of  Funding  Numbers 


Program    Element    Number    Project   No        Task  No        Work  Unit  Accession  No 


elude  Security  Classification)  Latent  Factor  Models  and  Analyses  for  operator  Response  Times 


i  Author(s)  Gaver,  D.  P.  and  O'Muircheartaigh,  I.  G 


)f  Report 

al 


13b  Time  Covered 
From  To 


14   Date  of  Report  (year,  month.day) 

1990,  September 


1 5  Page  Count 


nemary  Notation  The  views  expressed  in  this  paper  are  those  of  the  author  and  do  not  reflect  the  official 
position  of  the  Department  of  Defense  or  the  U.S.  Government. 


Codes 


JP 


Subgroup 


1  8    Subject  Terms  (continue  on  reverse  if  necessary  and  identify  by  block  number) 

Latent  factor  models;  extreme-value  distribution;  Weibull  distribution;  bootstrap; 
maximum  likelihood  estimation;  operator  response  times 


ct  (continue  on  reverse  if  necessary  and  identify  by  block  number 

models  are  presented  for  the  response  times  of  different  operators  to  different  tasks  where  response  is 
oy  one  or  more  cues  provided  by  the  system.  One  model  for  the  log-response  times  is  a  mixed  or  latent 
>del  with  unequal  case  fixed  effects  and  variances.  The  other  model  for  the  log-response  times  is  a  non- 
log-extreme-value  model.  Procedures  for  estimating  the  parameters  by  maximum  likelihood  are 
1.  The  models  are  used  to  analyze  response  time  data  from  simulator  experiments  involving  nuclear 
ant  operators  performing  certain  safety-related  tasks.  The  Findings  of  the  models  are  critiqued  and 
Dns  to  risk  analysis  are  sketched. 


bution/Availability  of  Abstract 
lassified/unlimited  same  as  report 


D 


DTIC  users 


2  1      Abstract  Security  Classification 

Unclassified 


:  of  Responsible  Individual 

iver 


22b    Telephone  (Include  Area  code) 

(408)  646-2605 


22c  Office  Symbol 

OR/Gv 


I  1473,  84  MAR 


83  APR  edition  may  be  used  until  exhausted 
All  other  editions  are  obsolete 

i 


security  classification  of  this  page 

Unclassified 


LATENT  FACTOR  MODELS  AND  ANALYSES 
FOR  OPERATOR  RESPONSE  TIMES 

D.  P.  Caver 

Department  of  Operations  Research 

Naval  Postgraduate  School 

Monterey,  CA  93943 

I.  G.  O'Muircheartaigh 

University  College 

Galway,  Ireland 


0.  INTRODUCTION,  BACKGROUND,  AND  SUMMARY 

There  are  many  situations  in  which  an  operator  (single  individual,  or  group 
or  crew)  is  confronted  with  a  somewhat  complex  task  that  must  be  accomplished 
within  prescribed  time  limits.  The  task  actually  often  initially  requires  diagnos- 
tic steps  followed  by  action.  In  some  cases  the  diagnostic  steps  are  stimulated 
by  a  cue  event,  leading  to  probing  actions  intended  to  reveal  the  correctness 
of  a  tentative  diagnosis,  followed  by  observation  and  interpretation  of  system 
response,  in  turn  followed  by  viewpoint  revision  and  further  action.  While  it 
is  intriguing  to  attempt  to  model  response  in  such  detailed  terms,  this  paper 
does  not  embark  on  that  enterprise.  Rather,  we  provide  and  analyze  models  for 
the  resulting  overall  response  time  of  different  operators  to  different  tasks  where 
response  is  initiated  by  one  or  more  cues  provided  by  the  system.  Two  factor- 
analytic  models  are  presented  along  with  likelihood  estimation  procedures.  The 
latter  are  then  employed  to  analyze  data  sets  from  typical  exercises  conducted 
at  simulators  used  for  training  nuclear  power  plant  operators;  their  identities 
are  kept  anonymous.  (The  findings  of  the  model  are  critiqued,  and  applications 
to  risk  analysis  are  sketched.) 
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It  is  believed  that  similar  models  will  be  useful  for  summarizing  the  behavior 
of  operators  or  crews  in  other  situations,  both  military  and  otherwise.  For 
example,  application  to  military  tank  driver  performance  is  envisioned. 


1.  A  MIXED  OR  LATENT  FACTOR  MODEL  WITH  THE 
UNEQUAL  CASE  FIXED  EFFECTS  AND  VARIANCES 
(LOG  N  MODEL). 

Consider  this  linear  model  of  mixed  type: 

Yik  =  (i  +  uk  +  ufi  +  eik         z  =  1, 2, ...  J 

k  =  1,2,...  A'.  () 

where  Yxk  =  InTjfc  with  T,^  being  the  time  for  crew  i  to  respond  to  situation  k; 
fi  and  uk  are  fixed  constants  (effects),  and  a>,  ~  IIDN(0,(t5),  e,jfc  ~  IIDN(0,ct£), 
are,  respectively,  the  latent  random  component  that  "individualizes"  case  (in- 
dividual, crew,  etc.)  i,  and  the  random  variation  displayed  by  any  individual 
on  situation  (task,  problem,  etc.)  A'.  It  is  assumed  that  each  case  occurs  in 
conjunction  with  each  situation  (e.g.  a  person  confronts  a  particular  problem) 
just  once  in  the  data  set  to  be  modelled.  In  practical  circumstances,  some  such 
individual  interactions  may  be  missing  for  reasons  unrelated  to  individuals  and 
situations,  a  problem  that  is  deferred  for  the  present;  see  Appendix  A  and  B. 

As  implied,  the  model  described  may  well  be  of  interest  when  data  pertaining 
to  human  performance  are  to  be  analyzed,  but  should  also  be  of  use  elsewhere. 
The  K  tasks  or  items  are  allowed  to  have  their  own  fixed  response  properties, 
described  by  (^fc,<T£);  this  pair  will  be  referred  to  as  a  task  signature.  The  usual 
mixed  ANOVA  model  formulation  assumes  a\  =  a2,  constant  for  all  k  (see 
Scheffe  (1959)),  as  is  reasonable  when  measurement  error  is  represented. 

Note  that  because  of  the  assumption  of  possibly  unequal  a^'s  afixed-o;,  model 
cannot  be  usefully  estimated  by  likelihood.  Consequently  the  above  random- 
effect  model  has  been  introduced,  and  fitted  to  data.  As  will  be  apparent,  it 
is  possible  to  estimate  the  posterior  density  for  U{  using  Bayes'  formula  in  the 
style  of  empirical  Bayes;  the  mean  of  the  resulting  Normal/Gaussian  density 
£[w,|data,  log  N  model],  is  available  as  an  estimate  of  the  ith  crew  effect  if  so 
desired. 

That  the  above  setup  is  a  latent  factor  model  has  been  remarked  to  us 
by  Professor  T.  W.  Anderson;   see  Anderson  (1988)  Chapter   14  for  relevant 


(overage  in  the  Normal/Gaussian  case.  Brillinger  and  Preisler  (1983)  survey 
various  other  latent  factor  data-analytical  studies  in  non-Gaussian  settings;  this 
includes  a  detailed  discussion  of  a  latent  factor  Poisson  model  for  counting  data. 
Brillinger's  paper  is  interesting  in  that  it  suggests  examining  goodness  of  model 
fit  by  "uniform  residuals,"  a  procedure  considered  in  our  study  as  well. 

Fitting  the  LOG  N  model  by  likelihood  requires  iterative  calculations;  the 
setup  is  described  in  the  next  section.  In  case  one  wishes  to  "robustify"  the 
formulation,  perhaps  by  introducing  more  outlier-prone  specifications  such  as  the 
Student  t  or  Tukey  density  for  a;,  then  more  numerical  effort,  or  approximation 
is  required.  Use  of  the  Laplace  approximation  together  with  Gauss-Hermite 
integration  may  well  turn  out  to  be  useful;  see  Gaver  and  O'Muircheartaigh 
(1987),  and  Gaver,  Jacobs  and  O'Muircheartaigh  (1990).  In  a  later  section  a 
totally  non-Normal/Gauss  model  for  operator  response  times  is  introduced  and 
fitted. 


2.  FITTING  THE  LOG-NORMAL  MODEL  BY  MAXIMUM 
LIKELIHOOD 

In  the  model  (1)  the  individualizing  case  effect,  a>,  for  case  i  is  viewed  as 
a  latent  or  unobserved  rv  whose  effect  on  the  Y{k  observable  is  indirect.  What 
is  the  probability  distribution  of  Y^^k  =  1,2,... A"  in  terms  of  the  unknown 
parameters?  Clearly  it  is  multivariate  normal  since  u>,  and  e,k  occur  as  a  sum; 
the  density  for  case  (crew)  i  is,  by  conditional  independence,  given  a>t  , 


K  e-|(2/,t-M-^-u>,)2/^ 
/£(£•;/*,£,<*>!•)     =][ 7= 

,_.  y/ZTT(7k 


k=l 


K 
-\1>2iVik-l*-i>k-Wi?tol  (2) 

e     *=i 


K 

(v^)An^- 

To  obtain  the  unconditional  density  of  y   remove  the  condition  on  u>,: 

The  calculation  needed  ("completion  of  the  square"  in  the  exponent)  can  be 
expeditiously  performed  as  follows.  Recognize  that  the  exponent  is  quadratic  in 

L>;  put 


A'i  being  independent  of  u>.  To  find  u>j,r2,  and  A',  differentiate  (2.3)  re  w  and 
equate  coefficients  of  u  and  1: 

(—  -      \        K 

^P-  =  £  (»«  -  /*  -  "*  -  «)  /***  (5) 

so,  from  the  cu-term, 

i/t*  =  y;i/*i  (6) 

fc=i 


wliil<-  from  t  he  1   tei  m 


ft 


u,/r2  =  ^2(yik  -  n-vk)l°l, 


(7) 


giving 


K 
w<  =  Y, (yik  -n-*>k)Wk  =  y,. -  {n  +  V.) 

k=\ 


(8) 


where  Wk  =  (l/oj)  r2  =  (1/ffJ)  /E/=i     l/^f?  %■  and  i/.  are  thus  W^-weighted 
averages.  Substitution  into  (4)  gives 


ft 


Kt  =Y/(ytk-»-vk-u>)  l°l 


k=i 


ft 


=  £ 

k=i 


Vik  -V-Vk-Yl  (Vik  ~  V  ~  Vk)  Wk 


k=\ 


l/°k 


K 


=  £[(**-*■) -to-"-)]  iM- 
fc=i 


Now  from  (2) 


Jo 


a  familiar  convolution,  from  which 


oo       c-i(57i-w)2/T2e-iw2/'2 


(9) 


rfu;,  (10) 


/k.U,;/^,^2,^)  = 


e-iA-,e-ip,)2/(r2+^) 


r-       /     A' 


(v^)    n^  v^  y^r^ 


^=i 


Thus  the  likelihood  of  n,v,a  ,  a2  is  proportional  to 


«=i 


n  **)  \A2 + ^ 


A:=i 


(ii) 


(12) 


and  hence  the  log-likelihood,  /,  may  be  expressed  as 
2 /(/*,£, £2,0£;^)  =  /    lnr2- 

K  I  I  (13) 

fc=i  t=i      t=i 


Differentiation  gives  these  estimates: 
9/  7 

— ;  =  °  =  £  [(^-l/.-)-K-^)](i/^),  (H) 


d{uk  -  i/.) 


1=1 


SO 

«=i  .=i  k=i  (15) 

=  (?.*-?..)• 

Next, 


dl  I  '      t 


0(/z  +  I/.)  r2  + 


(16) 


w    t=i 

SO 


*!+!/.  =  ?..  (17) 

These  of  course  closely  resemble  conventional  ANOVA  estimates. 

When  estimating  variances  it  is  convenient  to  reparameterize  in  terms  of 
precision:  pk  =  1/<tJ,    p  =  1/r2  =  £jfc=i  l/^jt     =  Ef=i  Pit- 
Then 


«  =  o  =  -L  +  J-  -  ^4  _  A?(t)  _  LiiHiZf „ 

0P*  P       P*       l/P  +  ^5  (l/p  +  a^f 


i  =  l 


and 


Next 


Aj(fc)  =  £[(»»-y,-.)-(y.*-p..)]2  (19) 

=  1 

/ 
<P)2  =  J2(Pi)2-  (2°) 


1=1 


91  =  o  =  -/^  +  7J2>U, 


a<rj  "  r2+<7j       (r'  +  ffj)2' 


yields 

i+,2  1^2  1     £     ^.f.  (22) 

^  t=l 

introduction  of  (21)  into  (18)  simplifies  the  latter  to 

-  =  7   A?   (*)  +  -,  *=l,2,...,fl\  (23) 

Pk       1  P 

The  system  (15),  (16),  (22)  and  (23)  must  be  solved  iteratively.  Begin  by  simply 

fitting  as  if  a\  =  a2  to  estimate  /i(l),^(l),u;(l),  and  obtain 

i  / 

^(1)=7ZT  E    (».-*  -  A(l)  -  Afc(l)  -  "»(1))2  (24) 

from  which  compute  ^(2)  =  (l/a£(l))  /  (e£i  i/^/2)-  Next  calculate  */fc-i/.(2) 
y.jt(2)  -  j/..(2)  using  W^(2)  in  (15),  and  fi+v(2)  =  y..(2)  from  (16).  It  is  now 
possible  to  evaluate  A/(fc;2)  from  (19),  and  (uJ)2(2)  from  (20),  and  hence  pk(2) 
and  p(2)  from  (23),  after  which  <r£(2)  from  (22).  Now  recompute  Vp)t(3)  = 
(1/(7^(2))  /  (2H/L1  l/^/2)  —  Pfc(2)/p(2),  and  so  repeat  iteratively  until  conver- 
gence is  achieved.  A  solution  procedure  based  on  Newton-Raphson  iteration 
has  also  been  obtained;  agreement  of  the  two  procedures  is  generally  good. 


3.  LOG-EXTREME- VALUE  MODEL  (THE  LOG  EV  MODEL) 

An  alternative  to  the  previous  model  that  may  be  attractive  is  the  following 
setup: 

(a)  T{k  is  distributed  according  to  a  two-parameter  Weibull;  then  it  follows 

mathematically  that 

(b)  Ylk  =  In  Tlk  has  the  extreme-value  distribution: 

1  -  exp  [-0,  exp  ((ytk  -  i)k)  /£*)],  with  probability  density  function 

fY,k(y,k;Vk,tk\0,)  = 

exp  [-0,  exp  ((ylk  -  T]k)  /&)]  0,  exp  ((ylk  -  r]k)/(k)  — 


(25) 


Note  the  occurrence  of  parameter  0,,  which  is  intended  to  represent  crew 
effect,  i.e.  0,  is  a  way  of  individualizing  crews  comparable  to  the  action  of  u,  in 
the  previous  model.  Values  of  0,  are  viewed  as  randomly  selected  latent  factors 
as  were  the  u>x  values.  The  nature  of  the  0,  contribution  differs  from  u>,  in  this 
model:  whereas  in  the  LOG  N  model  u,  acted  purely  additively  (on  the  log  scale) 
to  affect  the  center  (mean  of  logged  response  times)  in  a  manner  common  to  all 
tasks,  in  the  LOG  EV  model  it  can  be  seen  that  logged  times  are  represented  as 

1^  =  */*  +  &[(- In  0,-) +  *«•*],  (26) 

€ik  having  standardized  extreme  value  df.  For  the  present  model,  (25)  or  (26), 

E  [Yik\$i]  =  Vk-  0.5772^  -  6  In  0,  (27) 

Var  [Yik\9i]  =  ?-£  -  IMA^l  (28) 

6 

which  permit  initial  parameter  estimation  by  moments  and  facilitates  compari- 
son with  the  results  of  alternative  models.  Expression  (26)  implies  that  responses 
to  tasks  are  affected  differentially:  the  greater  the  natural  variation  in  perform- 
ing a  task  by  crews  (measured  by  £k  for  task  k)  the  greater  the  "average"  effect 


on  task  duration  due  to  crew  effect.  This  is  a  specific  form  of  interaction  be- 
tween  crew  and  task  effects  that  may  (or  may  not)  be  reasonable  in  particular 
circumstances. 

Conditional  on  0,,  the  crew  fs  response,  K,-,  on  K  different  tasks  has  joint 
density  function 

K 

fy.  (&;  a^^O  =  II^.»  (^  *»&'  ft),  (29) 

where  conditional  independence  is  assumed.  In  order  to  obtain  the  unconditional 
joint  density  of  response  Yj  remove  the  condition  on  dt  by  integrating  out;  this 
step  corresponds  to  (10).  Thus 

fy,  (&■;  2>i)  =  **.  {«"'■*'.■*} 4,  (3°) 

where 

K 

C^exp  ((**-%)/&)  (31) 

fc=i 

4  =  exp  (£  (ytt  -  %)/&)  f[  (1/fc)  (32) 

\fc=i  /  fc=i 

The  above  model  closely  resembles  one  introduced  by  Crowder  (1985)  and 

Crowder  and  Kimber  (1989).  However,  ours  deals  with  the  log  time,  and  hence 

is  a  location-scale  model  that  more  closely  compares  to  the  additive  log-normal 

model,  although  the  %  is  not  generally  a  mean,  nor  is  fa  a  standard  deviation. 


4.  GAMMA  VARIATION  FOR  9 

A  search  for  mathematical  tractability  suggests  that  variation  in  6  be  de- 
scribed by  a  gamma  density: 

„      r-O/0W)1//3~l     1      r      (I    l\  ,„* 

9~e        mm    ~0  =  Gam\0>v  (33) 

10 


and 


so  E[#]  =  1  and  Var[#]  =  /?,  and  from  which  the  joint  density  of  observations 
by  crew  i  is 

The  likelihood  associated  with  /  independent  crews  is 

_(T(K  +  l/P)aK\'+r(      1      \K  +  U0 


J  =  lnZ     =  IK\np  +  I\n(T(K  +  l/P)/T(l/P)) 


i  i  (36) 

-  (A'  +  1//3)£>(1  +  Pa)  +  £></,, 

t  =  l  :=1 

After  arrangement  and  re-parameterization  so  that  fa  =  ln£*.  the  log-likelihood 
becomes 

A'-l  /  IK  K 

I  =  /£  ln(fc  +  IIP)  -  (A  +  l/p)J^\n(l  +Pa)  +  J2   £(lK*  -  »»)<*P(-&)  "  7E^ 

fc=0  i=l  t  =  l    fc=i  fc=l 

(37) 
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5.  FITTING  THE  LOG-EXTREME  VALUE  (LOG-EV)  MODEL 
BY  MAXIMUM  LIKELIHOOD 

To  obtain  the  maximum  likelihood  estimates  of  the  parameters  we  iteratively 
solve  the  following  equations,  for  which  k  =  1,2,...,  A'  throughout: 

dl/dp  =  0, 

dl/drik  =  0,  (38) 

dl/dfa  =  0. 

One  Newton-Raphson  iteration  only  is  applied  to  each  equation,  after  which 
the  entire  process  is  repeated  until  convergence.  Typically,  only  two  or  three 
repetitions  are  required. 

We  record  the  derivatives  needed  for  the  above  process. 


ai  '  ! 

—     =     -(A'  +  l//?)^(ct/(l  +  /3cI))  +  (l//J2)^ln(l  +  /?ct) 
^  «=i  :=i 

K-\ 

-I   Y,  1/0(1  +  */*)  (39) 

k=0 

dl  ' 

5-     =     (A'/3+l)exp(-^)^er'V(l+/3ct)-/e-^  (40) 

°Vk  ,=1 


where  rlk  =  (yjk  -  %)/&,  a  residual;  finally 


dl  1  1 

=  (K0  +  1)  £  r,,er'V(l  +  fa)  -  £>,■*  ~  I-  (41] 


d<j>k 
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The  second  derivatives  are 
d2l 


dp 


2 


/  / 

=  IK/132  +  (K  +  1//3)  £  c'/(!  +  ^)2  +  (V02)  £  «.•/(!  +  ^) 
t=i  1=1 

-(2//?3)  x;  in(i + pa)  -  (i/p4)  x;  c* + upy2 

t=l  fc=0 


H 


-(2///33)£  (k  +  l/P)-\ 
k=o 


i 
■  2<f>k 


(42) 


=     (l  +  K(3)e-z^J2  {(er'V(l  +  ^,))[(/?er'V(l  +  /3ct))-l],    (43) 


t=i 


=     (!  +  #/?)£  {/?(r»«r«"*/(l  +  ^))a-er,'*»'.-*(l  +  rl-*)/(l  +  i9ci)} 


i=i 


+  !>.•*■  (44) 

!=i 

In  order  to  use  the  inverse  of  the  Fisher  information  matrix  to  provide  stan- 
dard errors  of  the  parameter  estimates  all  cross-partial  derivatives  are  required; 
we  omit  recording  these  in  the  interest  of  brevity;  the  expressions  may  be  ob- 
tained from  the  authors.  Our  numerical  experience  has  been  that  standard 
errors  obtained  from  Fisher  information  tend  to  be  too  small,  as  judged  from 
bootstrapping  approaches  next  to  be  described. 

6.  BOOTSTRAPPING 

A  modern  alternative  for  obtaining  standard  errors  and  approximate  confi- 
dence limits  is  the  parametric  bootstrap  of  Efron  (1979;  esp.  Remark  K,  p. 25). 
This  procedure  has  recently  been  applied  to  failure  data  in  the  context  of  the 
Challenger  disaster  by  Dalai,  Fowlkes  and  Hoadley  (1988)  and  goes  as  follows: 
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put0  =  (4>l,H,fi,i'2,"i'K;<>'i,--  -^A-)'n  Model  LOG  N,  and  (/5;  7ft,  . . .  ,r/A-;  0i,...,<fa) 
in  Model  LOG  EV.  Note  that  £  =  (9U...,0P)  here  denotes  a  generic  parameter; 
it   bears  no  direct  relation  to  the  ith  crew  effect  in  our  LOG  EV  model,  (25). 
Then  our  procedure  is  this: 

(a)  Estimate  0  from  data;  the  result  is  £(0),  the  point  estimate  of  the  param- 
eters. 

(b)  Provisionally  adopt  0(0)  as  the  true  value  in  the  parametric  model,  in  the 
present  case  (1)  or  (25). 

(c)  Simulate  B  independent  data  sets  (bootstrap  samples)  from  the  model 
evaluated  at  1(0):  {Yik(b),  i  =  1,...,/,  k  =  1,2,..., A';  6=  1,2,.  ..,£}. 

(d)  Compute  estimates  of  0  for  each  sample,  obtaining  the  bootstrap  estimates 
{0(b) J)  =  1,2,..  .,#}  =  {0(B)},  the  bootstrap  distribution  of  0. 

(e)  Present  relevant  statistical  summaries  of  marginal  and  joint  distributions 
of  {0(B)}:  e.g.  use  as  standard  error  of  0(0)  components  the  corresponding 
standard  deviations  of  the  bootstrap  estimate;  use  as  confidence  limits 
upper  and  lower  percent  points  of  the  bootstrap  sampling  distributions, 
suitably  adjusted.  We  present  numerical  illustrations  in  the  next  section. 

(f )  The  same  procedure  can  evaluate  standard  errors  of,  and  confidence  limits 
for  predictions  from  data:  in  the  present  case  prediction  of  the  probability 
that  a  response  time  exceeds  any  given  value  is  evaluated  in  terms  of  the 
model  evaluated  repeatedly  at  bootstrap  parameter  estimate  values;  see 
Dalai  et  al  (1988  )  for  an  example.  See  Section  9  for  an  example  in  the 
present  context. 

Bootstrapping  methods  suggest  themselves  for  comparing  the  adequacies  of 
different  models  for  fitting  and  predicting  from  specific  data  sets.  Specifically, 
bootstrapping  may  assist  in  choosing  between  two,  or  more,  candidate  models. 
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In  the  present  setting  one  may  wish  to  predict  the  probability  of  non-success, 
i.e.,  of  response  time  exceeding  some  time  window  of  duration  /,  P(t;  9).  Models 
A  and  B  (e.g.  our  LOG  N  and  LOG  EV  options)  are  estimated  obtaining  £4(0) 
and  0B(0).  Then  generate  bootstrap  samples  for  A  and  B  using  £4(0)  and 
9B(0)  respectively,  resampling  to  estimate  the  mean-squared  error  of  prediction 
when  Model  i  is  used  to  predict,  given  that  the  data  comes  from  Model  j;  here 
{(i,j)}  =  {i4,/4;  A,B;  B,B;  B,A}.  Prefer  the  model  whose  use  minimizes  the 
maximum  estimated  mean-squared  error  of  prediction.  An  alternative  strategy 
is  to  prefer  prediction  from  the  most  conservative  model:  the  one  predicting  the 
greatest  risk;  see  Section  9  and  Draper,  Hodges,  et  al  (1987). 

Another  option  for  residual  examination  is  to  compute  estimates  of  the  ex- 
pected log  response  times  associated  with  each  Task/Crew  combination.  Since 
crew  effects  are  random,  we  estimate  them  in  specific  cases  by  their  posterior 
means. 

For  the  LOG  N  model,  examination  of  (10)  reveals  that  the  posterior  density 
of  u>,  is  AT  [(ljJt2)  I  (1/r2  +  \/al)  ,  1/  (1/r2  +  1/V2 )]  .  We  substitute  in  the 
rule's  for  the  various  parameters  to  estimate  in  a  particular  case: 

tfi  =  (wi/(r)2)/(l/ra  +  l/a*) 

and  then  from  (1),  (8),  (15),  (16) 

Vik  =  (/i+iO  +  K1^.)  +  "i  =  y..  +  (y.k  -  y..)  +  [(</,  -  y..)l(ff\  /  [l/(f2)  +  l/al 

(45) 
where  all  averages  are  suitably  weighted.  Note  that  the  above  formula  for  the 
mean  acts,  in  effect  as  if  a  preliminary  hypothesis  test  for  homogeneity  of  crews 
is  being  applied:  if  cr2  is  very  small,  giving  evidence  that  all  crews  are  the  same, 
then  the  estimate  y^  ~  y.k,  the  (weighted)  task  mean  for  each  crew.  On  the 
other  hand  if  <r2  is  very  large  then  y^  ~  y.k  +  (?/:.  -  2/..)>  the  task  mean  modified 
by  the  estimated  effect  for  crew  i.  The  effectiveness  of  such  a  smooth  transition 
when  pooling  data  was  noted  by  Mosteller  (1947). 
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For  tlic  LOG  EV  model  t ak<-  the  expectation  of  (27)  with  respect  to  the 
(estimated)  posterior  density  of  0,,  which  is  Gammafe,  +  1/(3,  K  +  \/(i)  from 
(.'!.'{)  and  (34).   Using  the  first  two  terms  of  the  asymptotic  expansion  we  find 

V,k      =  E[yik) 

=  fjk  -  0.57726  -  £k  [in  (a'  +  1//?)  -  In  (c,  +  1/fl)  -  0.5/  [k  +  1/p) 

(46) 
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7.  EXAMPLE  DATA  ANALYSES 

The  previous  models  have  been  used  to  analyze  response  time  data  from  sim- 
ulator experiments  involving  operators  performing  certain  safety- related  tasks. 
In  Tables  1  and  2  appear  actual  data  from  two  such:  System  L  and  System  D. 
It  is  noted  that  certain  task-crew  combinations  are  missing,  some  because  of 
simulator  failure.  Two  procedures  were  adopted  for  dealing  with  these  cases: 
(1)  values  were  imputed  by  an  EM-like  process;  see  Little  and  Rubin  (1987); 
alternatively,  (2)  a  likelihood  approach  was  taken  that  simply  omits  such  values 
from  the  analysis  by  setting  to  unity  the  likelihood  contribution  associated  with 
a  cell  having  a  missing  response.  Both  approaches  can  be  useful;  the  former 
leans  more  heavily  on  model  correctness.  The  analyses  reported  here  emphasize 
the  use  of  a  simple  imputation  procedure;  the  incomplete  data  results  are  also 
reported  in  the  summary  tables. 

In  addition  to  "true"  missing  values  there  are  observations,  here  marked 
with  asterisks,  that  were  judged  to  result  from  operators  following  non-standard 
response  strategies.  It  was  judged  to  be  useful  to  analyze  the  data  both  with, 
and  without,  including  such  values;  when  omitted,  those  nonentries  in  the  data 
table  were  treated  as  missing  values,  entries  imputed  as  above  or  treated  as 
missing  values,  and  analyses  made  using  LOG  N  and  LOG  EV  models. 

Results  of  fitting,  along  with  bootstrapped  standard  errors,  appear  in  Tables 
3  and  4  for  System  L,  and  Tables  5  and  6  for  System  D.  Both  tables  exhibit 
main  effects  (/i  +  Vk,Vk)  and  scale  parameters  (o^and^*.)  for  LOG  N/EV  mod- 
els computed  under  A:  missing  values,  and  non-standard  strategy  values,  both 
treated  as  literally  missing,  using  methods  of  the  Appendix,  and  alternatively 
with  entries  imputed,  and  B:  only  the  "true"  missing  values  treated  as  missing. 
The  imputation  process  used  was  iteration  based  on  standard  two-way  ANOVA 
with  fixed  task  and  crew  effects.  This  is  a  convenient  crude  approximation  to  a 
proper  EM  algorithmic  approach,  Little  and  Rubin  (1987).  In  addition,  random 
crew  effect  variance  parameters,  <r£  and  (3  respectively  for  the  two  models,  were 

17 


estimated.  The  resampling-refitting  parametric  bootstrap  supplied  the  standard 
errors;  see  Efron  (1979)  and  Dalai,  Fowlkes,  and  Hoadley  (1988). 

It  is  noted  that  for  System  L  the  fixed  effects  (/z-f  Vk,Vk)  under  A  and  B  agree 
closely,  with  the  exception  of  Tasks  2  and  9.  Data  for  Task  9,  in  Table  1,  exhibits 
8  out  of  18  missing  values,  and  a  further  4  non-standard  strategy  values,  all  of 
which  are  far  in  excess  of  other  times  for  that  task.  Data  for  Task  2  also  exhibit 
substantially  many  missing  values  and  non-standard  times,  the  latter  having 
resulted  from  operators  following  non-standard  procedures  and  hence  yielding 
times  more  lengthy  than  the  other,  acceptable,  values.  The  noticeable  differences 
are,  however,  still  within  2  bootstrap  standard  errors.  The  corresponding  scale 
effects  (e.g.,  log  task  time  standard  deviations)  for  Tasks  2  and  9  behave  in 
corresponding  fashion,  increasing  by  factors  of  2  to  3  if  the  non-standard  times 
are  included. 

Similar  behavior  occurs  for  System  D,  although  here  the  exceptions  occur 
for  Tasks  4  and  8.  There  are  fewer  missing  and  non-standard  times  reported 
for  System  D  than  for  System  L.  Relatively  large  changes  occur  in  the  standard 
errors,  as  well  as  in  main  effect  and  scale  parameter  estimates,  when  several 
missing  or  non-standard  times  are  encountered,  and  these  are  treated  differently 
in  the  analysis. 

In  order  to  check  for  the  effect  of  imputation,  and  also  for  that  of  apparent 
correlations  between  certain  task  times  the  analyses  were  re-run  for  System  L 
omitting  Tasks  2,  4,  and  9.  The  results  appear  in  Tables  7  and  8.  Although 
specific  numerical  values  are  changed,  the  general  pattern  remains  quite  similar: 
the  new  numbers  are  quite  often  well  within  a  standard  error  of  the  estimates 
that  utilize  data  from  all  tasks. 
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8.  MODEL  CRITICISM  VIA  RESIDUAL  EXAMINATION 

In  order  to  examine  the  overall  fit  of  the  models  to  data  it  is  useful  to  con- 
duct some  form  of  residual  analysis.  We  have  chosen  first  to  judge  the  overall 
degree  of  fit  by  computing  and  summarizing  uniform  residuals  as  described  by 
Brillinger  and  Preisler  (1983).  In  general  for  this  procedure  one  estimates  model 
parameters  9,  from  data  and  then  examines  the  estimated  probability  integral 
transformation  of  the  data,  utilizing  the  fitted  model:  if  the  model  is  correct 
then  the  latter  should  closely  resemble  the  uniform  distribution.  In  Figures  1 
through  8  we  display  plots  and  summary  statistics  for  such  estimated  proba- 
bility integral  transforms  of  the  present  data  set.  We  also  exhibit  the  result 
of  bootstrapping  once:  each  model  was  allowed  to  create  one  set  of  bootstrap 
sample  data  utilizing  the  fitted  parameter  values;  these  values  were  then  treated 
like  raw  data,  and  were  then  probability  integral  transformed  and  the  results 
plotted  and  summarized. 

The  results  have  different  implications  for  the  appropriateness  of  the  models 
for  the  two  data  sets.  The  left  uniform  plot  of  Figure  1  shows  decided  non- 
uniformity  of  residuals  when  the  raw  data  is  fitted  by  LOG  N  using  the  methods 
of  the  Appendix;  however,  if  a  bootstrap  sample  is  generated  using  the  fitted 
model  parameters  the  results  are  far  more  uniform.  This  strongly  suggests  that 
the  basic  model  is  inappropriate.  A  similar  implication  is  obtained  by  examining 
Figure  3,  the  residuals  of  which  are  associated  with  LOG  EV.  Figures  2  and  4 
describe  the  residuals  when  imputation  of  missing  and/or  non-standard  values 
is  conducted.  Notice  that  uniformity  of  residuals  of  the  fit  of  the  raw  (plus 
imputed)  data  is  greatly  enhanced.  This  is  not  surprising  since  imputation  is 
based  on  presumption  of  model  correctness,  and  the  missing  and  non-standard 
values  are  imputed  using  the  presumed  model.  The  same  general  behavior  is 
observed  when  data  from  System  D  are  fitted  by  the  two  models  in  various  ways. 
Here,  however,  the  departure  of  the  residuals  from  uniformity,  as  shown  on  the 
left-most  residual  plots  of  Figures  5  and  7,  seems  less  pronounced  than  is  the 
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ca  ;e  for  System  L.  Again,  imputation  improves  the  uniformity  of  the  residuals 
and  the  apparent  fit.  We  conclude  that  the  log-additive  models  are  more  likely 
to  be  trustworthy  for  system  D  than  for  System  L. 
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9.  MODIFICATION  OF  THE  MODELS  FOR  INITIAL  DELAY 

In  the  nuclear  plant  simulator  exercises,  and  doubtless  for  other  applications 
as  well,  certain  task  response  times  may  begin  after  a  cue  different  from,  and 
later  than,  the  actual  initiating  event.  That  is,  the  operational  cue  that  triggers 
response  may  occur  some  time  after  a  possible  original  cueing  event  (e.g.,  the 
first  evidence  of  nuclear  plant  abnormality).  Unfortunately,  the  time  of  the  op- 
erational cue  is  not  usually  recorded  in  simulator  practice,  and  so  response  time 
data,  which  may  use  initiating  event  time  or  some  other  well-specified  event 
time  for  reference,  tends  to  exhibit  an  initial  delay.  Such  delays  may  be  inferred 
from  plant  and  initiating  event  information,  or  by  examining  the  actual  response 
time  data.  Note  that  ignoring  such  delays  if  they  are  appreciable,  i.e.,  fitting 
a  2-parameter  LOG  N  or  LOG-EV  model  when  a  3-parameter  specification  is 
more  appropriate  can  importantly  change  the  estimated  parameter  values,  par- 
ticularly the  LOG  N  o>  or  LOG  EV  £*.,  the  measures  of  within-task  variability. 
Specifically,  if  the  delay  is  ignored  and  our  current  models  employed  uncritically 
when  a  delay  7^  >  0  is  required  then  the  estimated  values,  &k  and  £k,  will  be 
biased  downwards  (under-estimated),  sometimes  quite  significantly. 

To  illustrate  the  effect  of  a  rough  accounting  for  delay  examine  the  Task  10 
data  for  system  L  in  Table  1.  The  minimum  value  is  1402  and  the  maximum 
is  3450,  so  there  is  an  appreciable  delay  associated  with  beginning  the  task 
as  compared  to  the  variability  of  the  response  times.  If  the  data  is  taken  at 
face  value  and  our  LOG  N  and  LOG  EV  models  fitted,  then  Table  3  exhibits 
£10  =  0.28  considerably  smaller  than  that  for  other  tasks.  If,  however,  a  rough 
adjustment  is  made  for  delay  by  subtracting  1000  from  each  Task  10  response 
time  data  value  then  the  simple  standard  deviation  estimate  for  logged  (times- 
1000)  for  Task  10  becomes  0.59,  far  more  similar  to  other  task  response  time 
standard  deviations. 

At  present  the  LOG  N  and  LOG  EV  programs  fit  2-parameter  models,  so 
accounting  for  delay,  7^,  must  be  done  off-line.    A  formal  approach  to  the  es- 
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timation  problem  for  the  Weibull  model  is  given  by  Smith  and  Naylor  (1987); 
that  paper  also  n>fW<  n<<     other  relevant  articles. 
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10.  RISK  CALCULATIONS 

An  important  application  of  the  LOG  N  and  LOG  EV  models  is  to  risk 
analysis:  it  is  desired  to  estimate  the  probability  that  a  task's  response  time 
occurs  within  a  particular  time  window,  for  in  that  case  (some  aspect  of)  the 
threat  has  been  averted.  We  will  refer  to  the  probability  that  response  time 
exceeds  the  time  window  as  the  human  interaction  risk  associated  with  the  task. 

It  is  easy  to  make  point  estimates  of  the  required  probabilities  using  both 
existing  models:  one  simply  replaces  the  model  parameters  associated  with  the 
task  of  interest  by  their  maximum  likelihood  estimates.  A  natural  way  in  which 
to  handle  crew  effect  is  simply  to  remove  the  crew  condition  u>  or  6  by  integrating 
it  out  with  respect  to  the  appropriate  Normal/Gauss  or  Gamma  estimated  prior. 
In  order  to  assess  the  effect  of  a  particular  crew  on  the  risk  it  is  necessary  to 
calculate  the  posterior  density  for  that  crew  and  then  integrate  out  on  u  or  9  with 
respect  to  that  posterior.  It  may  well  be  of  interest  to  compare  the  estimated 
risk  as  it  depends  upon  which  crew  is  in  place  when  an  initiating  event  occurs 
in  order  to  assess  the  effect  of  individual  crews  directly  on  risk.  This  is  not  done 
here. 
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11.   RISK  UNCERTAINTY 

In  order  to  assess  the  the  uncertainty  inherent  in  the  risk  estimation  one  may 
once  again  bootstrap.  The  procedure  is  this: 

(a)  Estimate  6  from  data  to  obtain  0O: 

LOGN:  <?0  =  (<^,/i»2) 
LOGEV:  fc  =  0M,f) 

(b)  Provisionally  adopt  0O  as  the  true  value  in  the  parametric  model. 

(c)  Simulate  B  independent  data  sets  from  the  model  with  parameter  0q  : 
{yik(b),  with  6=  1,2,. ..,B}. 

(d)  Estimate  $  from  each  sample  of  (c).  Obtain  {6(b),  b  =  1,2,. . . ,  B}\ 

(e)  Compute  P{T.k  >  wk\8(b)}  =  P{ln  T.k  >  In  wk\0(b)}; 
=  P{Yk>w'k\0(b)}  =  rk(b),  where 

wk  =  In  wk. 

The  risk  associated  with  Task  k  estimated  from  bootstrap  sample  b(b  = 
1,2, ... ,  B)  for  each  model,  that  is,  calculate  the  probability  that  the  time 
window  wk  is  exceeded  using  the  bth  set  of  parameters  estimated.  For  our 
two  models  and  b  =  1,2, .. .  ,B),  and  also  6  =  0,  the  original  estimate,  thus 
becomes,  respectively,  for  the  two  models, 

flt(6;LOGN)  =  l-#(^Jfiffi±M»V  (47) 

V  \l*im  +  *2(»)  / 

Rk(b;LOG  EV)=  (l  +  /3(6)exp  (w'k  -  f,k  (b))  /£(6)j  (48) 
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The  following  table  exhibits  risk  calculation  for  prescribed  windows;  the  lengths 
chosen  are  for  illustration  only.  Calculations  have  been  made  both  without  non- 
standard strategy  values  (A)  and  including  non-standard  strategy  values  (B); 
any  missing  or  left-out  values  were  imputed  as  before. 

Examination  of  the  point  estimates  in  Table  shows  that  there  is  considerable 
similarity  in  the  orders  of  magnitude  of  the  risks  given  by  LOG  N  and  LOG  EV. 
However,  LOG  EV  values  are  consistently  below  those  for  LOG  N,  as  are  the 
corresponding  confidence  limits.  It  is,  thus,  more  conservative  to  adopt  LOG  N 
model-based  risk  estimates  than  to  use  those  based  on  LOG  EV,  or  equivalently 
the  Weibull  model.  The  exception,  Task  9,  is  probably  traceable  to  the  many 
missing  values  exhibited. 

(f )  Present  statistical  summaries  of  marginal  and  joint  distribution  of  Rk(b): 
e.g.,  use  as  standard  error  for  the  original  risk  estimate,  Rk(0),  the  stan- 
dard deviation 


Sr"  =  \  B^rit(Rk{b)-RkY  (49) 

Note  that  the  overall  risk  associated  with  a  particular  initiating  event  de- 
pends upon  the  risks  associated  with  all  tasks  (human  interactions)  associated 
with  response  to  the  event.  These  risks  may  well  be  dependent  probabilities, 
at  minimum  because  all  tasks  presumably  confront  the  same  crew.  To  handle 
the  dependency  induced  by  a  common  crew  one  can  assume  conditional  inde- 
pendence, multiply  risks  conditional  on  crew  to  obtain  the  risks  associated  with 
joint  events,  and  then  integrate  out  with  respect  to  the  crew's  posterior  density. 
The  calculation  may  be  performed  explicitly  for  the  LOG  EV  model  since  a 
closed-form  elementary  function  expression  for  the  extreme- value  survivor  func- 
tion exists;  no  such  simple  calculation  can  be  carried  out  for  the  LOG  N,  but 
numerical  procedures  are  always  available. 

25 


Examination  of  the  results  suggests  considerable  similarity  between  the  risks 
calculated  using  the  LOG  N  and  LOG  EV  models.  It  is,  however,  noticeable 
that  for  all  tasks  except  Tasks  3  and  9,  LOG  N  predicts  a  slightly  greater  risk 
than  does  LOG  EV,  and  generally  with  slightly  larger  standard  error.  Of  course 
in  most  cases  shown  the  risks  associated  with  the  window  values  in  our  example 
are  too  high  to  be  realistic;  the  windows  were  chosen  for  illustration  only. 

The  above  calculations  have  been  carried  out  using  parameter  values  ob- 
tained under  imputation.  Thus  the  apparent  similarity  of  risks  across  models  is 
probably  overstated,  and,  in  view  of  the  suspicion  cast  on  the  fits  to  system  L 
data  (see  Section  8)  we  should  treat  the  System  L  risks  with  caution. 
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12.  SUMMARY  AND  CONCLUSIONS 

In  this  article  we  have  suggested  and  shown  how  to  fit,  by  maximum  likeli- 
hood, two  models  for  operator  response  times.  The  fits  of  the  models  to  two  sets 
of  (actual)  data  are  displayed  and  compared.  Uncertainties  are  assessed  by  para- 
metric bootstrapping.  Complications  involving  missing  values  and  consequent 
lack  of  balance  are  dealt  with  by  direct  likelihood  computation  as  well  as  by  a 
simple  form  of  imputation.  Finally,  the  fitted  models  are  applied  to  estimate 
the  risk  of  exceeding  (hypothetical)  time  windows;  associated  uncertainties,  i.e., 
standard  errors  and  confidence  limits,  are  obtained  by  bootstrapping. 

We  view  this  work  as  a  pilot  or  feasibility  study  intended  to  illustrate  and 
explore  possibly  useful  approaches  and  methodology  to  an  important  area.  Out- 
standing problems  remain:  the  models  put  forward  were  chosen  for  their  abilities 
to  account  for  some  aspects  of  the  real  situation  and  for  relative  tractability,  but 
many  other  forms  could  be  conjured  up,  fitted,  and  applied  to  infer  risk,  as  de- 
fined here.  It  is  somewhat  interesting  to  find  that  risks  estimated  from  the  same 
data  using  the  different  models  agree  rather  closely;  it  is  not  unlikely  that  the 
agreement  will  suffer  if  the  windows  are  increased  so  as  to  achieve  much  smaller, 
and  presumably  more  realistic  risk  values. 

The  bootstrap  standard  errors  and  confidence  limits  warn  that  although 
the  present  models  tend  to  agree  in  their  risk  assessments  the  uncertainty  is 
still  rather  large,  even  if  wrong-model  or  structural  uncertainty  is  ignored;  see 
Draper,  Hodges  et  al  ( 1987)  for  discussion.  In  order  to  reduce  the  uncertainty  of 
estimation,  e.g.  to  reduce  standard  error  size,  it  is  often  proposed  to  aggregate  or 
pool  data,  either  from  similar  tasks  in  the  same  environment  (plant),  or  for  the 
same  task  across  "similar"  environments.  Both  procedures  are  worthy  of  investi- 
gation, but.  will  be  credible  only  if  suitable  adjustments  are  made  to  reduce  bias. 
Adjustments  can  be  carried  out  by  using  models  resembling  the  types  suggested 
here,  possibly  enhanced  to  include  regression  terms  "explaining"  responses  in 
terms  of  measured  and  qualitative  crew  and  plant  characteristics  ("performance 
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shaping  factors"  is  the  jargon  in  certain  risk  assessment  circles).  To  date,  for- 
iii. J  adjustment  attempts  by  regression  have  been  inconclusive  but  are  a  form  of 
insurance  that  should  be  included,  and  validated  to  the  greatest  extent  possible, 
if  aggregation  or  pooling  is  contemplated,  particularly  across  plants.  Inter-plant 
variability  may  be  appreciable  because  of  variations  in  management  philosophy 
and  style. 

In  general,  it  seems  advisable  to  utilize  bootstrapping  as  extensively  as  possi- 
ble to  build  an  appreciation  for  the  variabilities  and  uncertainties  involved  when 
using  models.  Bootstrapping  that  uses  direct  re-sampling  as  originally  discussed, 
Efron  (1979),  seems  difficult  or  impossible  for  situations  such  as  are  described 
here  unless  vastly  more  data  becomes  available  and  better,  more  scientifically 
based,  models  and  true  replications  can  be  employed.  Consequently,  use  of 
parametric  bootstrapping  becomes  necessary.  However,  parametric  bootstrap- 
ping that  consumes  data  from  more  realistic,  and  elaborate,  models  and  their 
fits  to  simpler  structures  can  be  useful  and  informative.  But  such  exercises  can- 
not directly  substitute  for  data  obtained  under  truly  operational  circumstances, 
which  even  the  best  simulator  data  can  not  aspire  to  be. 
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APPENDIX  A.  MAXIMUM  LIKELIHOOD  ESTIMATION  FOR  LOG 
N  MODEL,  WITH  MISSING  VALUES 

Let  dlk  -  1  if  the  observation  at  (i,k)  [i.e.,  for  crew  i,  task  k]  is  available; 
otherwise  d,k  —  0.  Then  the  component  of  likelihood  associated  with  (Crew)  i 

is 


4fe,/^,-.)  =  n{-  -jjj 


>.KOk 


K 

e      *=> 


(V2*) 

Now  write 


(v«)*"n£.i»i" 


*'  1  /.-.  .x2 


«.*  (&*  -H-vk-u)    •  -j  =       —J +  A',. 

After  differentiation  re  u>, 


fc=i  a*  r^ 


T  (ytfc  -  n-vk-u)  = 


which  implies,  identifying  terms  of  order  1  and  u>,  and  writing 

A' 

E  4*  (!/•*  -  A*  -  uk)Pk  =  w,-/r? 
ik=i 

A 

E  «tap*  =  1/t? 

A 

Hence  r?  =  1/  £  dlkPk 
k=\ 

K  K  K 

"i  =  Ti    E  dik  (ytk  -n-  uk)pk  =  yL  -  p.  -  r f  £  dlkvkpk,  where  yt.  =  r f  £  dikyikpk 

and 

A 
A',  =  ^  rflJk  (yrt  -  fl-  „k  -  u.  f  pk , 
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Thus  re-write  the  likelihood  as 

>  rH(*-)2/^}e-K/^ 


J  y,  (lU'>li'&Z*t£)  =e  *Ki  J 


(X 


not*}  yfif+*u 


k= 
Hence  the  likelihood  assumes  the  form 


Taking  logarithms,  we  obtain 

'  ,-.2 


i  =  2\ul  =  -J2  *  -  E  ^rh>  +  Elnr*  -  E  n*ln<7*  -  Eln  {T>  +  <£) . 

1=1  1=1    «  +  aw      1=1  k=l  t=l 


=  -E*-i:^3+D»'f+i;«»i«»-i:i«('?+«i) 

1=1  1=1    «  "*"    w      1=1  Jt=l  1=1 

where 

«fc  =  E  rf«* 
i=i 

Differentiating  first  re   /*,  we  obtain 


^77  ~  ~~  ^      /)„ 


dKj  _      dux  _     1 
Hence 


M_         ^2u>,-(-l)=() 


V - =  0 


31 


Vi.  -  V  -  Tf  £  dik"kPk 

k=l 


{rf  +  al) 


=  0 


/  _j        /    \Vi.  -rf^dikp^kj 


fife.  -  TfY.dxk^kPk\ 


**  =  ?        W  +  «» 


/E(*  +  «2) 


-l 


Next,  differentiate  re   i//: 


I       '        w 


Y  d*k  (y>k  -p-vk-  Uif  pk 


-  2  Y  dik  (Vik  -  P~  »k  ~^x)Vk\  -hk  ~ 
k 


dui 


Hence 


m_ 


=  -rfdupi 


Y  \  2  Y dlk  (y*k  -  V-  uk  -Ui)pk  \-6\k  +  rfdupij  +  2   2    '    2    -rfdtlpi    >  =  0 


1=1     v         k 


Thus, 


/ 


Y  du  \  TidUPl  Y  dtk  (yik  ~  V-Vk-  W,)pjfc  -  dxk  (yit  -  fJ,  -  Vi  -  LJt)pi 2  '   '   2dupi  >  = 


t  =  l  V  k 

Further  simplification  results  in 


_  Y  dtl  {  Ti  Y  d'k  (y'k  ~  V  ~  vk  ~  U{)pk  -  (yu  -  p  -  i/i  -  Ui) 


U>iTf 


rU°l 


=  0 
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Hence 


v\  ^2  dtl  =  Yl  dt< 


i-\ 


~Ti    J2  d^  (&'*  -P-Vk-  "i)Pk  +  (Vxl  ~  fi  ~  Wj)  +   ,    2',     \. 

ft  \Tt  +  ai) 


if  ?_ 

=> vi -  —  Y2 du \~ Ti  Yl dik (&* - /*  - "* - Q*)pk  +  (yw - /* - wj)  +  ,  2r',u,'n 

n<  ,=1     I       ft  (Tl  +  °l) 

We  now  differentiate  re  pi,  noting 


dpi 
du)t 


£  dikPk 
k 


-du  =  -rfdu 


-z-j-  =  rfdu  (yil  -  fi  -  i/i)  +  1^2  dik  (ylk  -  ft  -  uk)pk  (-r^d,/) 


dp/ 


=  T-?d,-j  <  ftj  -  fi  -  vi  -  t?  Y^  dlk  (yik  -  fi  -  uk) pk  \ 


=  rfdti(y,i  -  \i-vi  -u>i) 


OK, 

dpi 


=  du  (yu  -  \i-vi-  u>,)2+2^  dlk  (ylk  -  /z  -  vk  -  ujt)pk  {-T?du  (yit  -  p,  -  v\  -  cj,)} 


=  du  (yik  -  n  -  vk  -  u>i)  <  (yu  -  \i  -  v\  -  w,-)  -  2r,2  J^  dlk  (yik  -  /x  -  uk  -  Wj ) pft  f 


Noting  also  that 


r2  J^  rf,ft  ( l/.it  -P-Vk-  Wj )  Pfc  =  r2  J^  d,-ft  ( j/,-fc  -  n-  Vk)Pk  -Wi  r2  X!  ^'*P*  =  w'  ~  w«  =  ° 

ft  ft  ft 


it  follows  that 


Bhi  2 

dp, 
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di 


'Mr? 


=  Y]  dxl  (yu  -  n-vi-  u;,)2  -      2   '  '       (yt,  -  p.  -  u,  -  Qt) 

dpi     fri  (r, +  ai) 

(-l)w,2  (1  -  rfdu)        1    /      4j  x    ,   n,       \(-T?d«) 


tf  +  ol)' 


+  ~2    \  ~Ti  dill   + 7-2— oT 


Hence, 


dpi 


J^dii  \  -{Vxi  -  H-  "l  -Wi) 


t=i 


2w,T2 

T  +  ^2. 


u;2r4 
1    t 


(tf  +  O* 


rf  + 


rf  +  ai 


Pi 


Pi  =  ntj  <  Y2d'1 


1=1 


(</«/  -  A*  -  "/  -  Wj )  I  y,-/  -  ^i  -  1//  -  o>j  + 


2^,r2 


u>2r4 
t    1 


^ — 1-  r"  - 

(r,2 + -If 


t'  +  °\ 


Finally,  difTerentiating  re  <72 ,  we  find 


1L  =  T_A y      '      =0 


which  implies 


'  /".2  T2  '  1 

The  following  iterative  procedure  was  employed: 

(a)  Use  imputed  values  to  obtain  p(0)  ,  <72  (0)  ,7*  (0),p/t  (0)[A:  =  1,...,A"] 

(b)  Solve  for  p(l)  using  (A-l)  Updated, 

(c)  Solve  for  0t  (1)  using  (A-2)  Update  u>,- 

(d)  Solve  for  pi  using  (A-3)  Update  u>,-,  r2 

(e)  Solve  for  <r2  using  (A-4) 

Repeat  (b)  through  (e)  until  convergence. 
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